Drift Laws for Spiral Waves on Curved Anisotropic Surfaces 
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Rotating spiral waves organize spatial patterns in chemical, physical and biological excitable 
systems. Factors affecting their dynamics such as spatiotemporal drift are of great interest for par- 
ticular applications. Here, we propose a quantitative description for spiral wave dynamics on curved 
surfaces which shows that for a wide class of systems, including the BZ reaction and anisotropic 
cardiac tissue, the Ricci curvature scalar of the surface is the main determinant of spiral wave drift. 
The theory provides explicit equations for spiral wave drift direction, drift velocity and the period 
of rotation. Depending on the parameters, the drift can be directed to the regions of either maximal 
or minimal Ricci scalar curvature, which was verified by direct numerical simulations. 
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Introduction. Spiral waves of excitation have been ob- 
served in diverse chemical, biological and physical sys- 
tems [THS] . They organize spatial patterns of excitation 
and underly important processes such as morphogenesis 
of a social amoeba [7] , some forms of neurological dis- 
ease [8] and cardiac arrhythmias [3|9]. In many cases, 
the dynamics of spiral waves is of great interest because 
it determines the overall behavior of the system. One of 
the most important aspects of the dynamics is the spa- 
tiotemporal drift of spiral waves. The drift of spirals can 
determine the type of cardiac arrhythmia [S]; it has also 
been observed in the BZ reaction [TOUT], CO oxidation 
on a Pt surface [3] and biological morphogenesis [7j . Cur- 
rently, several sources of spiral drift have been identified, 
including tissue heterogeneity [T2"Hr5] . an external elec- 
trical field [TMTB] . spatially varying anisotropy [T§H2"3"] 
and surface curvature [21] ■ The latter is highly relevant 
as most of real excitable media have complex geometries 
which may include curved domain boundaries, e.g. the 
walls of human atria have a complex curved shape. 

Essential questions regarding drift of spiral waves on 
a surface are: What determines the drift direction and 
velocity of the drift? How is spiral wave drift affected by 
anisotropy of the medium? Can the parameters of drift 
be predicted from general properties of 2D spiral waves? 
Some of these questions were addressed previously in the 
kinematic approach 23 - 25J . However, the kinematic ap- 
proach used there is valid only for spirals with a large 
core, i.e. where front-tail interactions are absent. How- 
ever, curved surfaces with anisotropy were never studied 
before, despite their usefulness for cardiac applications. 

Here, we propose a theory of spiral wave drift on 
curved anisotropic surfaces based on a gradient expansion 
around the spiral wave solution. We derive equations for 
drift of a spiral wave on a surface of arbitrary shape with 
anisotropy. We show that the drift velocity is given by 
the gradient of the so-called Ricci curvature scalar (RCS) 
of the surface, which arises as the generalization of the 
Gaussian curvature of a surface. The coefficients in our 
equation for drift velocity are explicitly obtained from 
the properties of the two-dimensional spiral wave solu- 
tion in an isotropic planar medium using response func- 



tions [26 29 . Interestingly, depending on parameters the 
drift can be directed to regions of the surface with either 
lowest or highest RCS. As the spatial distribution of the 
RCS can be easily computed (see Appendix [lj , the pro- 
posed theory can predict the regions which will attract 
or repel spiral waves in each particular situation. Note 
that for anisotropic diffusion, the extrema of RCS do not 
necessarily coincide with the places of extremal surface 
curvature. An example is provided in Fig. [T] 

We verify our theory by a direct comparison with nu- 
merical simulations and show that the derived equations 
predict with high accuracy the trajectories of spiral wave 
drift on curved surfaces with significant anisotropy. 

Reaction- diffusion (RD) equation. We start from the 
RD equation in terms of Cartesian coordinates x 1 . 
Anisotropy is built in through the diffusion tensor , 
whose eigenvalues are proportional to the squared con- 
duction velocities along the local material axes: 

d t u(r,t) = (D^i^djPui^t)) + F(u(f,t)). (1) 

This PDE describes how a state vector u of the system 
changes due to local processes F(u) and anisotropic dif- 
fusion. The constant, dimensionless matrix P allows to 
exclude some state variables from diffusion. We now de- 




FIG. 1: Curved anisotropic surfaces colored according to their 
Ricci curvature scalar, (a) Paraboloid z = — 0.05(a; 2 + y 2 ) 
with isotropic diffusion, (b) Same paraboloid with anisotropic 
diffusion (Dl/Dt = 4) with fiber angle a = (n/40)(x + y). 
White bars indicate the local fiber direction. 
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rive the analogue of Eq. ([TJ) on curved surfaces with 
isotropic or anisotropic diffusion. 

Isotropic diffusion on a surface. Any smooth sur- 
face can be parameterized as x l (s A ), (i <G {1,2, 3}, A £ 
{1,2}), where the s A form a curvilinear coordinate sys- 
tem. The gradient and divergence operators in the dif- 
fusion term of Eq. ([!]) should therefore be expressed us- 
ing the metric tensor Gab that is induced on the sur- 
face by the coordinate transform x l — > s A , with compo- 
nents Gab = dAX l dijdsx^ . In case of isotropic diffusion 
= D S lj ), RD systems are thus described by [23] 
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GG AB d B Pu) +F(u) 
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Anisotropic diffusion on a surface. When diffusion is 
anisotropic, the diffusive current for a given diffusion ten- 
sor D equals J = —D ■ gradPu . Transformation to 
surface coordinates brings J A = — D AB dsPu, where 
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{1,2,3}) (3) 
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Taking the divergence will give -^Oa ( 

as the diffusion term for curved, anisotropic surfaces. 
From Eq. ph, however, it follows that det(D AB ) = 
det(GAs))" 1 det( J D y ). Here, det(Z>^) is the product of 
the diffusivities along the local material axes and in most 
of the cases assumed to be constant [TS]. Thus, following 
work by Wellner et al. [5U] as in [3T] [35] we can write for 
curved surfaces with anisotropic diffusion that 



dtu = -^_d A (ygg AB d B Pu) + F(u), 



(4) 



with the metric tensor still given by the matrix inverse 
of the diffusion tensor, albeit in surface coordinates s A : 



9ab = DoiD-^AB = D {D-%d A x l d Bi 



(5) 



The constant factor Do has been included to make gAB 
dimensionless; we also put P = DqP. 

Drift equations. First, we introduce Riemann normal 
coordinates around the spiral's rotation center at t = 
|33j . Such coordinates explicitly reveal how curvature 
affects the metric in the region close to the origin, in 
terms of the Riemann curvature tensor Rabcd- This 
tensor contains second order spatial derivatives of the 
metric 33J ; its trace is the Ricci curvature scalar 1Z men- 
tioned above. Furthermore, we associate an order A to 
each spatial differentiation of the metric tensor, as we 
are working in the regime of slowly varying anisotropy 
and small Gaussian curvature K q of the surface in com- 
parison with the spiral's core size. Hence, the fact that 
Riemann normal coordinates are locally Euclidean can 
be written g EF = Sef + 0(A 2 ). Therefore, the exact so- 
lution u(p A ,t) can be approximated by an unperturbed 
spiral solution u (p A ), i.e. u = u + u, where u = 0(A 2 ). 
Hereto, the solution is made to rotate at a frequency 
w around the spiral's center at position X A . In Ap- 
pendix [2j we show in detail how to bring this ansatz 



in the RD Eq. Q expressed in Riemann normal coordi- 
nates. With L = PA + ujode + F'(uo), we finally obtain 

Lu u + X a 8au q + (u) — wo)<9eu = S. (6) 

From the full derivation, S = TZS R - d A TZS A R + C(A 4 ) is 
found to contain spatial derivatives of the RCS and the 
unperturbed spiral solution Uq: 



u — Prd r u 



(7) 



S^ R = -^PpArdrUo + ^PpAd^Uo + ^P^dAU^ 
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The components of spiral drift can be found by project- 
ing Eq. (p3]) onto the critical adjoint eigenmodes Y - of L 
[TT1 12"51 [291 |M] , which are known as response functions 
[2"7] and were observed to be strongly localized around 
the spiral wave's tip. This property ensures that dy- 
namics of the full spiral wave can be essentially captured 
by expansion in Riemann normal coordinates close to its 
tip. Our procedure first delivers the instantaneous laws 
w - w = (Y e | S), X A = (Y A | S), where the bracket 
notation denotes an inner product in function space, i.e. 
(f | g) = J R2 dSi H g. If the drift and spiral core radius 
are small compared to the distance over which curvature 
and anisotropy change, one may average over one rota- 
tion to find the net spiral wave drift as 



dt<f> = LOo + qaU + OiX 4 ), (8a) 
d t X — — qigrad7\!. — q^n x gradTvL. (8b) 

where n is a unit normal vector to the surface, the grad 
operator taken with respect to the metric ^ and 
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(9a) 
(9b) 



The Eqs. @, are the main analytical result of 

this Letter. 

Discussion. The law ( 8b ) for spatial drift is strikingly 



similar to the law of motion of electrons in a solid mate- 
rial [35], where d t X = —figi&dcj) under an electric field 
E = —grad <j>. We will thus henceforth call qi the 'spiral 
mobility'. For positive mobility qi, the spiral wave will 
descend the gradient, ending up in a locus of minimal 
RCS, while for negative qi the spiral will drift to the re- 
gion with maximal RCS value. Spiral waves, however, 
exhibit also a second component of drift with coefficient 
#2 which makes them drift under an angle tan -1 qilq\ 
with the direction of the gradient of RCS. However, qi 
cannot influence whether the spiral drifts to higher or 
lower RCS. If the sense of spiral rotation is reversed, q\ 
remains the same, while qi switches sign. Since the RCS- 
induced drift contains third order spatial derivatives of 
the diffusion tensor, it is different from the metric drift in 
[15] . While the proportionality constants for the metric 
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drift equal the filament tension [53], we have found no 
simpler expression for the spiral mobility gi, g2- 

Let us now consider how the laws ([8]) apply to an 
isotropic diffusion system, as for example the BZ reac- 
tion. Here, K/2 = R 12 12 = K a , so the RCS is simply 
twice the Gaussian curvature of the surface. Eqs 
then confirm results obtained by Zykov et al. 
which were obtained in the kinematic approach. Our 
theory demonstrates that those results hold not only in 
the large core regime but for any stationarily rotating 
spiral wave. Also, it was suggested in [53] that in the 
equal diffusion case (P = I) the drift component parallel 
to the gradient of 1Z disappears, i.e. q\ = 0. However, 
in our theory this is not the case, and q\ can have any 
value depending on properties of the spiral wave. Our 
prediction is numerically confirmed in appendix [4| 

For cardiac tissue and other anisotropic RD systems, 
Eq. ( |8b[ ) is the first analytical expression that captures 
the dynamics of spiral waves on a curved anisotropic 
surface. Spiral drift is shown to be related to the gra- 
dient of the RCS which depends both on curvature of 
the surface and tissue anisotropy. When the local di- 
rection of maximal diffusivity is known in the medium, 
we may consider a smooth anisotropic surface with con- 
stant principal diffusivities Dl = AlDq^Dt — (It Do, 
where Dl > -Dt > 0. The respective eigenvectors of the 
diffusion tensor will be denoted et, e~r ; in the context 
of cardiac tissue, et is known as the local fiber direc- 
tion. In a local Euclidean frame one may then write that 




FIG. 2: Drift coefficients 90,91,92 induced by Ricci scalar 
curvature in Barkley's model 36 (b = 0.19, e = 0.025, P = 
diag(l,0)) for varying a. Signs of 90,92 for counterclockwise 
rotation. 



In accordance with Eq. (8a), this explicitly shows that 



-,AB 



spiral waves rotate faster on sphere-like surfaces, as can 
be expected from the angular deficit, thereby extending 
the results of [25] to non-uniformly curved surfaces with 
anisotropic diffusion. Note, that since the rotational cor- 
rection is small compared to ujq, only the spatial drift is 
further discussed here. 

To study both positive and negative mobility, we 
took a = 1.1 or a = 1.3, for which Eqs. fcty respec- 
tively predict (91,92) = (—0.102,2.652) and (91,92) = 
(0.643, 0.357) if the spirals rotate counterclockwise. In a 
first numerical experiment, we considered an anisotropic 
plane with linear fiber rotation as in |38| . i.e. et = 
d T S AB + (dt - d T )e£ ef . Hence, with a unit cosae^ + sinae^ with fiber angle a(x,y) = B (x + y). 

For such anisotropy, a direct analytical calculation gives 
1Z — 4(dt — cIt)B 2 sin 2a., which is color-coded in Fig. [3] 
The configuration does not possess isolated maxima or 
minima of the RCS: the local extrema are located along 
lines at — 7r/4 to the x axis. The minima occur at the 
fiber angle a = — 7r/4, while the maxima are found where 
a = 7r/4. We observe in Fig. [3^, that for positive mobility 
91 the spiral wave drifts towards minimal value of RCS, 
as predicted by our theory. In addition, we see a good 
correspondence of the real computed trajectory (green) 
and the trajectory obtained from integration of Eq.([8]) 
(black). For the negative mobility 91 in Fig. [3j3, we ob- 
serve only a small drift component towards the maximal 



normal to the surface, the drift law (8b I can be written, 



with V the gradient operator on an isotropic surface, 

X = - qi [d L ei (el ■ VTZ) + d T e T (e T ■ VK)] 

-q 2 Vd L d T eN x VK + <D(X 5 ). (10) 

For strongly anisotropic tissue (Dl S> Dt), the drift 
will thus appear to occur almost along the local fiber 
direction, as long as |gi| » |ga|< 

Comparison with numerical simulations. We verified 
the laws of motion Eqs. Q by direct numerical simu- 
lations. We used Barkley's reaction kinetics [35], where 
u = [u, v] T , F = [/ = e _1 u(l — u)(u— (v + b)/a), u — v] T 
and P = diag(l,A,). First we took e = 0.025, & = 
0.19, D v = and computed for varying model param- 
eter a the actual values of 90 , 91 , 92 by evaluation of 
the integrals Q using DXSPIRAL software [37]. Fig. [2] 
shows dependency of the coefficients 90, 91, 92 as function 
of the parameter a which determines the excitability of 
the medium (the higher values of a correspond to the 
higher excitability of the medium) . We see that the coef- 
ficient 91 is positive for most values of a, indicating drift 
to the lower values of RCS. However, in a medium with 
low excitability, the spiral mobility q± can be also neg- 
ative, making the spiral waves drift into the regions of 
higher RCS. For counterclockwise rotating spirals, the 
coefficient 92 is always positive and slightly decreases 
with a, while the coefficient 90 is negative and increases. 
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FIG. 3: Attractive and repulsive sites for spiral waves in the 
anisotropic plane of Fig. [Tjj, i.e. fiber angle a — (ir/4Q)(x + y) 
and Dl/Dt = 4. (a) Drift trajectories for Barkley's model 
as in Fig. [5] with parameter a — 1.3. (b) Same for a = 1.1. 
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FIG. 4: Spiral wave drift on a paraboloid surface; white bars 
indicate fiber direction. Drift trajectories for isotropic diffu- 
sion (yellow) and anisotropic diffusion as in Fig. [T]d (red) to- 
gether with theoretical predictions (black). Barkley's model 
as in Fig. [2] and parameters a — 1.3, A = 0.05, B = 7r/40 
(panel a) or a = 1.1, A = 0.025, B = tt/80 (panel b). 

value of RCS, as for this parameter value \qx\ « \q 2 \. 
Here too, however, the theoretical (black) and computed 
(green) trajectories almost coincide. 

In a second numerical experiment, the planar surface 
was replaced by the paraboloid z = —A(x 2 + y 2 ) with 



the same anisotropic properties as in Fig. [3] Its Gaus- 
sian curvature and RCS were shown in Figs. [1^,-b. We 
see that the observed drift trajectory (Fig. |4j red lines) 
is in close agreement with the theoretical predictions 
(black) for both cases of positive and negative q\. We 
also provide trajectories for spiral wave drift in absence 
of anisotropy, i.e. only due to Gaussian curvature of the 
surface (yellow lines). We see that for isotropic case the 
spiral indeed drifts away from the top of the paraboloid 
(qi > 0, Fig. [|a,), or slowly towards it (qi < , Fig. [IJd). 

Conclusions. We have developed an asymptotic the- 
ory that predicts the drift of spiral waves on general 
curved anisotropic surfaces. This drift is caused by a 
gradient of the Ricci curvature scalar and determined 
the coefficients relating this gradient and drift velocity 
using perturbation theory. The analytical results were 
quantitatively confirmed by numerical simulations. 
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Supplementary material 

1. Finding the metric and RCS for a surface 
with given shape and projected fiber angle 

We consider surfaces z = f(x, y) in the domain (x, y) € 
[-L/2,L/2] x [-L/2.L/2]. The surface is thought to 
contain fibers in the direction eX tangent to the surface, 
whose fiber angle is defined by tana = (el ■ ej)/(el • e%). 
For a prescribed angle a(x,y), one therefore finds that 
e~L = N [cos aelc+SYa.ae y + (cos a. d x f + sm.a d y f)e z ]; the 
factor N(x,y) is chosen to give el unit length. Now, we 
assume for a moment that the three-dimensional space 
is filled with copies of such surface in the direction of 
e" z . If diffusion along the local fiber direction occurs with 
diffusion coefficient Dl = Do^l, while transverse diffu- 
sion has Dt = DqcIt, the three-dimensional anisotropy 
is determined by the tensor 



D lj = D T 5 i] + (D L - D T )eie 3 L . 



(11) 



In the curved-space approach, a metric with contravari- 
ant components g 1 ^ = D^/Do is found in the three- 
dimensional space, with inverse gij. In our simula- 
tions, we chose to let the surface parameterization s A , 



A = 1,2 be 



x, s = y. From the transformation law 



9ab = d A x l gijdsx 3 , one finds 

ffii = 9xx + 2g xz d x f + g zz (d x f) 2 , 

922 = 9yy + fyyzdyf + 9 ZZ (dyf) 2 , (12) 

.9i2 = 9xy + 9xzd y f + g yz d x f + g zz d x fd y f. 

Thereafter, the metric components g AB are found as 
the matrix inverse of (g AB ). From the coefficients gAB, 
g AB , it is straightforward to compute the RCS using the 
Christoffel symbols T^ c [33] : 



r A _ 9 

1 BC — ~ 



AD 



(0b9cd + d C 9BD - d D g B c) , 



(13) 



K = g BC (d A Ti c - d c T 



A 

AB 



r A r D r A r D \ 
1 AD 1 BC 1 CD L AB) 



In index notation, our law of motion (8b) becomes 



X A = 



-Q.19 



-AB 



d B n-q 2 g- 1 / 2 e BA d B n + o(\ 5 ). (u) 



This expression is particularly useful for practical calcu- 
lations or numerical implementation. 

For the surfaces in Figs. [TJd and|4j numerical differen- 
tiation of g^B was used to obtain the RCS, with space 
step dx = 0.05. 



considered. The explicit expansion for the metric tensor 
in these coordinates can be found in differential geometry 
textbooks (e.g. [33]), yielding 



g EF (p\p 2 ) = S EF + ^Reabf(0,0)p a p b ~ 



(15a) 



12 



[d c R E ABF(0, 0) + 8cRfabe(0, 0)] + 0(X 4 ) 



In this expression, Riemann curvature components and 
its derivatives were evaluated at the center of rotation of 
the spiral wave solution, where p 1 = p 2 = 0. To find the 
metric components with upper indices, a matrix inversion 
is performed: 



9 EF (p\p 2 ) = 5 



EF 



1 



-R l 



A „B 



AB 



(0,0)p A p 



1 

12 



(15b) 



[8 C R E AB F (0, 0) + d c R F AB E (0, 0)] p A p B P c + G(A 4 ). 

Since we are dealing with two spatial dimensions only, 
we can use the identity 

Rabcd = (R-/2)(g Ac g BD - 9ad9bc) (16) 

whence, omitting terms of C(A 4 ), 

72.(0,0) 



Reabf (0,0) 



dcR 



C-n-EABF 



(0,0) 



2 

11(0, 0) 



d c n(o,o) 



($EB$AF — SefSab) 
-tEAtBF (17a) 

SefSab) 



(<>eb8af 



d c n(0,0) 



-tEA^BF 



(17b) 



Hereafter, we write the diffusion term of Eq. Q as 
^d A {^gg AB d B Vn) (18) 



19 



dA^i) g AB d B Vn + d A (g AB d B Vn 



Using the property g 1 l 2 d A g 1 ^ 2 



9 

\g EF d A 9EF, the 
right-hand side of Eq. ( 18 ) is found equal to d A d A Pu — 
■K(0,0)S R - a A 7e(0,0)S™ + C(A 4 ), with S R and S^ R 
given by Eq. Herein, dg = e A B p A d B , rd r = p d A 

and r 2 = S AB p^p B . 

Furthermore, we let the origin of the Riemann normal 
coordinates move along with the spiral wave's rotation 
center and rotate at the yet unknown rotation frequency 
lu, as in [31] . This requirement imposes a gauge condition 
on the field corrections u: 



2. Details of the drift equation derivation 

An essential part of our derivation is to express the 
diffusion term of the reaction-diffusion equation Q in 
Riemann normal coordinates. In these coordinates, the 
radial lines from the origin are geodesies of the space 



0. 



(Y x I u) 



0. 



(Y y | u) = 0. (19) 



Finally, the differentiation with respect to time r in the 
moving frame will generate convection terms |31j 



d T u = d t u — Lodgn — d t X d A u. 
which also appear in Eq. (J6|. 



(20) 
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Details of numerical simulations 



Evaluation of the coefficients qo,qi,q2 
using response functions 



The numerical values for the coefficients (i — 0, f , 2) 
displayed in Fig. [2] and listed as q\ h in the main text 
were acquired using an extension of the publicly avail- 
able dxspiral software in the following way. First, 
dxspiral was used to generate a standard spiral solution 
u for Barkley's model (a = 1.1, 6 = 0.19, e = 0.025, 
P = diag(l,0)) on a disc of radius R = 12.0, using a 
polar grid with N r = 240 and N g = 128. Thereafter, 
response functions Y , Y 1 , Y y were computed by the 
dxlin.c routine, whose details are given in [29] , The crit- 
ical eigenmodes and response functions were obtained in 
a complex basis, i.e. 



V 



(0) 



-d e u , V (1 



(1) 



(d x u Q - idyU ) , (21) 



w (o) = _ Y e , W (1) = - (Y x — iY v ) . 



Next, the overlap integrals (|9j were evaluated using the 
trapezoid rule, with terms S , S^ R given by Eq. (j7|). 

The coefficients shown in Fig. [2] for model parameter 
a > 1.1 were found in steps of 0.025 up to a — 1.4 by 
calculating new solutions u using the solution for the 
previous a as an initial guess. For each value of a, the re- 
sponse functions and the overlap integrals were evaluated 
as above in order to find qo, q\,q2 and 71,72- 



b. Forward simulation of the reaction- diffusion equation 
on curved anisotropic surfaces 

To check the validity and limitations of the theory, 
forward evolution of spiral waves was studied on curved 
anisotropic surfaces. Hereto, the reaction-diffusion equa- 
tion ([2]) was discretized using the finite difference tech- 
nique. With the purpose of studying generic surfaces 
whose shape is prescribed by z = f(x, y) in Cartesian 
coordinates, the curvilinear coordinates on the surface 
were taken to be s 1 = x, s 2 — y. That is, the function 
u^(x,y) would provide a top view on the field of the j- 
th variable of the spiral wave. A rectangular grid with 
dx = dy was taken. For the examples considered, we had 
d x f(0, 0) = d y f(0, 0) = 0, such that the finest spatial grid 
on the surface was obtained in the origin. This value also 
determined the largest time step allowed in our explicit 
Euler scheme; we chose dt = 0.9dx 2 /(4max(D L ,D T )). 

The diffusion term in Eq. |2| was discretized using 
a nine-point scheme, with the metric g the inverse of 
§ab from Eq. ( 12 ). For a given time t and state variable 



label j, we took 

^pd A (^-gg AB d B ui{x,y,t)) a - ]_ 

V9 dx 2 y/g{x,y) 

C m . n (x, y)u J (x + mdx, y + ndx, t). (22) 

m,ne{-l,0,l} 

Simple finite differencing yields the coefficients 
C m ,n(x,y), which were only computed at the start 
of the simulation and then stored. With h AB = ^[gg AB , 
they are 



Co,o(x,y)=- ^2 h11 \ x + m -^-^y ) 

m=±l V 1 ' 



(23a) 



h 22 I x, y + n 



n=±l 



dx 



C m fi{x,y) =h [x + m—,y 



dx 



E 

n=±l 



mn 



h 12 



x,y + n 



dy 



Co,n(x,y) =h 22 [x,y + 



dx 



(23b) 
(for to = ±1) 

(23c) 



E mn ,,, / dy \ 

— h [x + m—,y\, (for 77 = ±1) 



m=±l 

77777 , 12 ' dx 



C m , n (x,y)=—h l£,y + ?7 2 



mn, , 12 I dx 



(23d) 

(for 777, 77 = ±1). 



An overview of simulation parameters and grid size and 
resolution is presented in Tab. [Tj Before each simulation, 
a spiral wave was first created in a planar domain of 
larger size, same resolution and constant anisotropy equal 
to gAB at (x,y) = (0,0). The midpoint of the circular 
tip trajectory was determined, such that the standard 
spiral wave solution could be copied and centered on a 
suitable position in the anisotropic curved surface. This 
method allowed to reduce the duration of the transient 
regime and the associated drift, and therefore brought 
more control of the initial spiral wave position. 



Fig. a b A 


B L dx D v 


D L 


3a 1.3 0.19 


tt/40 30 0.1 





4 


3b 1.1 0.19 


tt/40 30 0.1 





4 


4a (red) 1.3 0.19 0.5 


tt/40 40 0.1 





4 


4a (yellow) 1.3 0.19 0.5 


40 0.1 





1 


4b (red) 1.1 0.19 0.025 


7r/80 80 0.1 





4 


4b (yellow) 1.1 0.19 0.025 


80 0.1 





1 


5 0.7 0.19 0.1 


40 0.1 


1 


1 



TABLE I: Overview of simulation parameters. All simulations 
had e = 0.025,A) = 1, A, = 1 and D T = 1 
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c. Numerical integration of the tip trajectory 

The dynamical laws Q predict the trajectory of the 
spiral's rotation center on curved anisotropic surfaces. 
The trajectory for a given initial position and surface 
shape and anisotropy can be found numerically by solv- 
ing the ordinary differential equation ( 8b I by going to 



index notation (14) first. This approach was taken for 



the trajectories on the curved anisotropic surface of Fig. 
[4j with time step dt — 0.05, values qi, q2 as predicted by 
dxspiral and g AB taken in coordinates s 1 = x, s 2 = y, 
i.e. given by the matrix inverse of Eq. ( fl2] ). 

The trajectories for the paraboloid with isotropic dif- 
fusion and anisotropic plane of Fig. [4j [3] were found as 
analytical solutions to Eq. ( 14 ) . 



4. Supplementary results 

Fig. [3] displays the drift trajectory of a spiral wave 
on a paraboloid surface with equal diffusion. Barkley's 
model was used for reaction kinetics, with a = 0.7, 
b = 0.19, e = 0.025 and P = diag(l,l). Although the 
kinematic approach in |24j states that q\ should vanish, 



it is clearly seen that the spiral wave drifts away from 
the top, in accordance with Eqs.(8b)-([9]), which yield 
(9i 5 ?2) = (0.855,-0.386) for a counterclockwise spiral. 



4 



Spiral tip 
Spiral center 
Spiral center (theory) I 




FIG. 5: Drift trajectory of a spiral wave on a paraboloid 
of revolution z = 0.1(a; 2 + y 2 ) in the equal diffusion case, 
showing a nonzero drift component along the gradient of the 
RCS. Colors indicate the RCS. 



